From XML dataΒΆ
Example on how to load your data if it is currently in xml format
import json
import dolfin
from fenics_plotly import plot
import pulse
try:
from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
from dolfin import Mesh, Function, Constant, DirichletBC, interpolate
# Mesh
mesh = Mesh(pulse.utils.mpi_comm_world(), "data/mesh.xml")
# Marker functions
facet_function = dolfin.MeshFunction("size_t", mesh, "data/facet_function.xml")
marker_functions = pulse.MarkerFunctions(ffun=facet_function)
# Markers
with open("data/markers.json", "r") as f:
markers = json.load(f)
# Fiber
fiber_element = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=4,
quad_scheme="default",
)
fiber_space = dolfin.FunctionSpace(mesh, fiber_element)
fiber = Function(fiber_space, "data/fiber.xml")
microstructure = pulse.Microstructure(f0=fiber)
# Create the geometry
geometry = pulse.HeartGeometry(
mesh=mesh,
markers=markers,
marker_functions=marker_functions,
microstructure=microstructure,
)
activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
activation.assign(Constant(0.0))
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
activation=activation,
parameters=matparams,
f0=geometry.f0,
)
# LV Pressure
lvp = Constant(0.0)
lv_marker = markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [pulse.RobinBC(value=Constant(base_spring), marker=markers["BASE"][0])]
# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
bc = DirichletBC(V.sub(0), Constant(0.0), geometry.ffun, markers["BASE"][0])
return bc
dirichlet_bc = [fix_basal_plane]
# You can also use a built in function for this
# from functools import partial
# dirichlet_bc = partial(pulse.mechanicsproblem.dirichlet_fix_base_directional,
# ffun=geometry.ffun,
# marker=geometry.markers["BASE"][0])
# Collect boundary conditions
bcs = pulse.BoundaryConditions(
dirichlet=dirichlet_bc,
neumann=neumann_bc,
robin=robin_bc,
)
# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
problem.solve()
(2, True)
# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, 0.2))
2021-11-23 16:52:20,054 [626] INFO pulse.iterate: Iterating to target control (f_26)...
2021-11-23 16:52:20,055 [626] INFO pulse.iterate: Current control: f_26 = 0.000,0.000
2021-11-23 16:52:20,056 [626] INFO pulse.iterate: Target: 1.000,0.200
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 78),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 252),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 389),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 450),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 511),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 572),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 641),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 762),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 823),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 892),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1005),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1110),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1275),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1328),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1405),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1519),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1589)],
[(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 74),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 76)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 248),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 250)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 385),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 387)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 446),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 448)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 507),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 509)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 568),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 570)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 637),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 639)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 758),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 760)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 819),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 821)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 888),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 890)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1001),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1003)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1106),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1108)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1271),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1273)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1324),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1326)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1401),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1403)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1515),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1517)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1585),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1587))])
# Get the solution
u, p = problem.state.split(deepcopy=True)
volume = geometry.cavity_volume(u=u)
print(f"Cavity volume = {volume}")
Cavity volume = 1.1409709245419595
# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, show=False))
fig.show()